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ABSTRACT 



Current numerical barotropic and baroclinic forecasting models have 
shown an inability to correctly represent the behavior of the planetary- 
scale waves. A mathematical model, based on the thermodynamic energy 
equation, has been proposed by Wiin-Nielsen as a means of predicting the 
behavior of the planetary-scale waves. The model is programmed for the 
Control Data Corporation 1604 digital computer to utilize the operational 
output of machine-processed data and analyses produced by the U. S. Navy 
Fleet Numerical Weather Facility (FNWF) , Monterey, California. The 
present form of the* program yields a 1-hr planetary-scale forecast at 
1000 mb, 850 mb, 700 mb, 500 mb and 300 mb in ten minutes. 
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1 . 



Introduction. 



Current numerical barotropic and baroclinic forecasts of the height 
of the 500-mb surface show an inability to correctly represent the 
behavior of planetary-scale waves. Early models caused these long waves 
to retrogress at speeds approaching the Rossby speed. This phenomenon 
is not observed in the atmosphere and much work has been done in an 
effort to correct this difficulty. Wolff [7] has shown that forecasts 
by the non-divergent barotropic model are more nearly correct if the 
planetary-scale waves are held stationary. Another approach to correct 
this defect has been to introduce a Helmholtz term into the vorticity 
equation. This is equivalent to adding a small amount of divergence. 
Wiin-Nielsen [5] pointed out that an analysis of baroclinic instability 
and phase speeds of waves in a baroclinic atmosphere leads to the con- 
clusion that the phase speed for large wave lengths approaches the 
Rossby speed in a non-divergent atmosphere. 

Burger [l] has shown, from scale considerations, that the vorticity 
equation loses its prognostic value and assumes a quasi-stationary 
character for waves whose lengths are of the same order of magnitude as 
the radius of the earth. Adopting this as a working hypothesis, and 
defining waves of this magnitude as planetary-scale waves, Wiin-Nielsen 
[6] developed a prediction model for the planetary-scale waves. 

This investigation is an attempt to provide a means of correctly 
predicting the planetary-scale motion based on the Wiin-Nielsen model. 
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2 . 



The Model. 



The vorticity equation for large-scale atmospheric flow in isobaric 
coordinates is 



a? 



V-V71* 






V} Vt 






(1) 



Burger [1] showed that this equation may be greatly simplified by order- 
of -magnitude considerations when applied to the planetary scale. 

Equation (1), after Burger, reduces to the stationary form 



yS\T ^ ^ V ‘ V “ O - (2) 

Wiin-Nielsen [6] demonstrated that equation (2) leads to instability on 
all scales, with large amplification rates for the smaller scale motions 
which may be present. The shorter waves were stabilized by adding the 
advection of relative vorticity into equation (2), leading to 



v-w? - £ v-v ■ o 

or (3) 

The coriolis parameter in equation (2) is modified to ^ in equation (3) 
to avoid mean -vorticity generation at any level. The inclusion of rela- 
tive vorticity leads to the delineation of a region of instability. In 
the meteorological range of vertical wind shear and wave length, the 
shortest e-folding time is of the order of 5 days. The e-folding time 
is the time required for the amplitude to increase by a factor of 0 

22 fctt 

Since the amplification rate is given by t- , the e-folding time 
is simply 

t e - \-/zVa . 
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The thermodynamic energy equation is derived in a straightforward 
manner from the first law of thermodynamics, applied to adiabatic flow s 



giving 



iQfVV-VQfV *[% UJ-^1 



• o. 



(4) 



This equation may be transformed to equation (5) by approximating 
by finite differences. 



& y 1 - o 



(5) 



where 



C r = 



R rvr BJ\ 

r W r c p ' 



( 6 ) 



To apply equations (3) and (5), the atmosphere is divided into 
9 layers by isobaric levels numbered 0 to 9, where the levels 0 and 9 
correspond to 200 mb and 1000 mb, respectively (see figure 1) . 
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Figure 1. Level arrangement for the model. 
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The vertical velocity is computed at even levels by applying the 



vorticity equation, (3), at the odd levels as follows: 



A = A-i + % X: V(f k „ <■ 0. (?) 

Here k is even and P is the pressure difference between even levels • 
The vertical velocity is assumed to be zero at 200 mb. 

The prognostic equations are obtained by applying the thermodynamic 
energy equation, (5), at the even levels, namely, 



s 



-X, -Vk 



'g; 



‘i 



Equation (6) may also be expressed as 



(8) 



BIT * SU-ll) (9) 

f. \ 

or 

s^l a ~ c ) 2 k 

= 3 V c r i • (10) 

Equation (9) is used to evaluate the stability parameter in a standard 
atmosphere at the even levels. These values are included in the program 
as constants. 

The lower boundary condition is established by noting that 



ui = = || * V- Vr ~/ p y / • (11) 

The second term of the right hand side of equation (11) is zero in 
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geostrophic flow. Since W~ at the surface is given by 



W= ^ .VWl , 

equation (11) may be rearranged to obtain equation (12a) 

-sr »r * 

or, approximately, 



r • V Ht 



1 



(12a) 



(12b) 



**t " ~ iT 'bp 

where the surface wind, % , has been approximated by . The 

difficulties embodied in equations (12a) and (12b) are rather serious 
since u) at the lower boundary is usually obtained by neglecting 



in equation (12a). Thus, it is likely that 



is a small differ- 
ence of two relatively large terms. In spite of these misgivings, the 
initial program was designed in accordance with Wiin-Nielsen 0 s model as 
approximated by equation (12b). 

The computations carried out during this investigation were made 
with a Control Data Corporation 1604 (CDC 1604) digital computer, with a 
core storage capacity of 32,768 words of 48 bits each. The input data 
were obtained from the bi-daily analysis of the height of the 300-mb, 
500-mb, 700-mb, 850-mb and 1000-mb surfaces by the U. S. Navy Fleet 
Numerical Weather Facility (FNWF) , Monterey, California. 

In order to utilize available data and computer subroutines from 
FNWF, the equations were converted to finite-difference form and scaled 
for fixed-point fractional arithmetic. 

A 1977-point, square-net, octagonal grid, circumscribed by latitude 
ION on a polar stereographic projection of the Northern Hemisphere, was 
employed to represent the data fields. The net spacing is 381 km at 60N 



where the projection is true. 
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The scaled, finite-difference equations corresponding to equations 
(8), (7), (10) and (12) are equations (13), (14), (15) and (16), 
respectively . 




The scaling conventions used in these equations are 



u) = 


A o' 7 

u>. Z 


(mb - sec 


? - 


w 


(sec ) 




o*2' 7 


(cm) 


P * 


? *2. 


(mb) 


m ~ 


in > 2' 


(non-dimensional) 


n* 


71-2 


(sec ) 


7 = 


/s 

z-2. 


(cm) 


cr a 


A ~ 11 

(T • 2 


2 -2 -2 
(cm sec mb ) 



(13) 



(14) 



(15) 



(16) 



(17) 
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(cm) 



(cm) 



Equation (6) represents the standard stability plus an anomaly due 



to the departure from the standard atmosphere. The quantity OT in 
equation (15) refers to the stability parameter at level k in a standard 
atmosphere. The remainder of the right hand side of equation (15) is the 
finite-difference form of equation (10), which was used to approximate 
the stability anomaly at each grid point. The point values of stability 
are averaged over each pressure surface to maintain consistency with 
respect to the energetics of the model. 

The vertical derivatives in the calculation of the stability anomaly 
are approximated in a conventional manner [3, pp. 47-49]. The height 
fields at the even levels (figure 1) are approximated by a simple linear 
relationship. It is assumed that the ratio of the height of a level to 
its adjacent level in a standard atmosphere is the same as the ratio of 
the heights of these levels in the real atmosphere. A test of the 
validity of this assumption indicates that the standard deviation of the 
actual ratios with respect to the standard ratios is two orders of magni- 
tude less than the ratios themselves. The approximation leads to a maxi- 
mum error of 9% of the height in the 54 cases examined, with an average 
error of less than 1%. The approximation is less accurate in the lower 
portions of the atmosphere than it is in the upper regions below the 
tropopause. 

Additional boundary conditions are imposed on the finite-difference 
approximation to the model. These include: 

a) setting the forecast height and thickness changes at the grid 



A 



boundary points to zero; 
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b) setting the relative vorticity and vertical velocity at the 
grid boundary points to zero. 

The order of computation during a time step is as follows: 

a) the previously prepared planetary-scale height analyses at the 
five odd levels (see figure 1) are entered into the computer 
and the height values at the five even levels are approximated; 

b) the stability parameter is computed at each grid point on 
levels 2, 4, 6 and 8 and is pressure averaged for each level; 

c) the vertical velocity is computed at each even level from the 
advection of absolute vorticity at the next lower odd level 
after establishing the vertical velocity at level 0 as zero; 

d) the thickness fields between adjacent odd levels are obtained 
from the input height fields and the advection of thickness is 
calculated; 

e) the advection of thickness and the vertical velocity term at 
each even level, except level 0, is combined to form the local 
thickness change; 

f) the height change at level 9 is determined; 

g) the calculated height and thickness changes are combined with 
the initial height at level 9 and the thicknesses of the upper 
layers to form the forecast height fields at levels 1 9 3, 5 9 

7 and 9. 

Forward-difference approximations are used in the initial time step. 

All subsequent time steps employ central differences. 
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3. 



Investigative Procedures and Results. 



An initial attempt at evaluating the model was made by providing 
the space-mean height fields at 300 mb, 500 mb, 700 mb, 850 mb and 
1000 mb as input. These fields are obtained by passing the analyzed 
height fields through four scans of the smoothing equation, 



In this instance, k = 1/8. This approach proved unsuccessful. An 
analysis by Holl and Clark [2] of the wave -component residuals shows 
that, after smoothing, 507** of the energy of wave number 16 remains in 
the space-mean pattern at latitude 45N„ This is sufficient justifica- 
tion to remove the FNWF space-mean, in its present form, as a suitable 
input for the model under consideration. 

The 1000-mb height-change forecasts produced by the approximations 
of equations (12a) and (12b) were unsatisfactory. The forecast changes 
had no pattern and were an order of magnitude larger than observed 
values. This result casts further doubt on the utility of equations 
(12a) and (12b). 

A second approach was to provide the model with input fields which 
were prepared by a hemispheric Fourier analysis of the height fields. 

The zonal mean and the components due to waves 1, 2 and 3 were extracted 
and added together. These summed fields provided a planetary-scale 
input to the model. Figure 2 is a sample of one of the fields used in 
the limited evaluation of the model. 

The first attempts with the Fourier input were not successful. 
Short-wave components were generated and first appeared in the mapped 
height fields after a 2-hr forecast. These components continued to grow 




( 18 ) 



in amplitude and progressed around latitude circles, as shown in figure 3 
for an 8-hr forecast. These short-wave phenomena are tentatively 
attributed to computational instability,, The time interval used in 
these unstable forecasts was one hour . 

The 1000-mb height-change forecast was not unstable; however, the 
earlier difficulties persisted. On this basis, the utility of equations 
(12a) and (12b) was rejected. In the interest of evaluating the main 
features of the model * a simple approximation to the lower boundary 
conditions was adopted. The 12-hr height changes at the 1000-mb surface 
were determined from the result of the Fourier analysis of the initial 
height fields. It was then assumed that the appropriate fraction of this 
change would be applied to the 1000-mb surface at each time step* Such 
an approximation leads to minor inconsistencies during the iterative 
procedure, but leads to a perfect 12-hr forecast at 1000 mb. 

The problem of possible computational instability has been 
approached on the basis of the general considerations outlined by 
Thompson [4]. The 1-hr time interval was halved, while retaining the 
approximation for the lower boundary. 

At the time of writing, 26 time steps have been calculated with no 
sign of short-wave phenomena, which appears to support the earlier sug- 
gestion of computational instability associated with 1-hr time steps. 
Figure 4 is a 12-hr forecast for the 500-mb level from 12Z 19 Novem- 
ber 1962. The errors associated with this forecast are shown in 
figure 3. The forecast 12-hr height-change field is figure 6 and the 
actual 12-hr height-change field is figure 7 . The analyzed 500-mb field 
for 00Z 20 November 1962 is figure 8. The pillow and root-mean-square 
error from figure 6 are very similar to those obtained from figure 7 
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(Pillow, P, and root -mean-square error, RMSE, are defined in Appendix I ) 
In figure 6, P = -2 feet and RMSE s 71 feet, while figure 7 gives 
P s 6 feet and RMSE s 74 feet, An examination of figure 5 shows a region 
of large errors in the vicinity of the North Pole (at the center of the 
figure) . This is to be expected since the features of the input data are 
not of planetary scale at these latitudes , The major errors, with a 
maximum of -480 feet, are associated with the slightly deepening low 
center in the northwest region of figure 2, The model forecasts a 
lowering of central height to 16,560 feeto In actuality, the height 
lowered to 16,770 feet at the position one grid point to the northeast 
Most of the other errors in the lower latitude regions are less than 
100 feet 0 These errors appear small enough to warrant further evaluation 
of the model. Appendix I is a tabulation of the pillow and RMSE of error 
fields generated when using persistence as a forecast on selected days of 
November 1962. These indicate the forecast improvement required of the 
Wiin-Nielsen model. 
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4. 



Suggestions for Future Study • 

The limited results described in section 3 are not yet adequate to 



evaluate the usefulness of the long -wave model 0 There remain several 
tasks to be accomplished before this may be established* 

a) Revise and streamline the present computer program to reduce 
the amount of time required to make the calculations for one 
time step. At present, this time is 5 minutes; however » this 
may be reduced, by more sophisticated tape handlings to about 
3 minutes. Further reduction might be made by packing data in 
the computer memory to minimize or eliminate the necessity of 
data transfer to and from magnetic tape, 

b) Conduct an extensive evaluation of the model on many sets of 
data to prove or disprove the adequacy of equations (3) and 
(5). 

c) If the evaluations of b) prove the model to be successful 9 a 
suitable lower boundary condition must be devised and tested 
in order to obtain operational planetary -scale forecasts. A 
further desirable feature for an operational forecasting model 
would be the elimination of the required Fourier analysis of 
the input data. The most likely possibility here may be the 
use of one of the new smoothing techniques presently undergoing 
test at FNWFo 
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Figure 2 - 500-mb height analysis of the Northern 

Hemisphere after Fourier analysis. 

12Z 19 November 1962. 
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Figure 3 - 500-mb planetary-scale height forecast 

for the Northern Hemisphere showing 
features evaluated a a computational 
instability. 8-hr forecast from 12Z 
19 November 1962. 1-hr increments. 
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Figure 4 - 12-hr 500-mb planetary-scale height forecast 

for the Northern Hemisphere from 12Z 19 
November 1962. 30-min increments. 
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Figure 5 - 12-hr 500-mb planetary-scale height -fore cast 

errors for the Northern Hemisphere from 
12Z 19 November 1962. 30-min increments. 
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Figure 6 — 12-hr 500-mb planetary-scale height-change 

forecast for the Northern Hemisphere from 
12Z 19 November 1962. 30-min increments. 
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Figure 7 



- TIP— hr 500-mb planetary-scale height-change 
analysis of the Northern Hemisphere from 
122 19 November 1962. 
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Figure 8 - 500-mb height analysis of the Northern 

Hemisphere after Fourier analysis* 

00Z 20 November 1962. 
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Figure 9 



Polar stereographic projection 
of the Northern Hemisphere to 
the scale of figs. 2 through 8 
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APPENDIX I 



VERIFICATION OF PERSISTENCE FORECASTS 
FOR SELECTED DAYS OF NOVEMBER 1962 



Pillow is defined as 






p = 






(I-D 



The root-mean-square error (RMSE) is determined as follows: 



a 

• ( 1 - 2 ) 

In equations (1-1) and (1-2) , ^ and represent the point 

values of the fields being compared . In this example, when preparing 
the data for these tables, the A-field was the analysis at the initial 
time and the B-field was at a subsequent time,, 



RMSE = 



x.-n* |W7 



£Ra-b)-p] 



Table I 
Table II 
Table III 



Index of Tables 

Pillow and RMSE for 12-hr forecasts 
Pillow and RMSE for 24-hr forecasts 
Pillow and RMSE for 36, 48, 60 and 72-hr forecasts 



Table format - Pillow / RMSE (values in feet) 
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Table I. 


Pillow and RMSE of persistence verifications 
for a 12-hr forecast during November 1962, 


at 5 levels 


Initial 


Time 






Level 










300mb 


500mb 


700mb 


850mb 


lOOOmb 


12Z 19 




22/99 


6/74 


0/57 


-2/50 


0/61 


OOZ 20 




-23/97 


-9/68 


3/54 


2/57 


-12/69 


12Z 20 




-11/88 


-3/71 


-3/63 


0/66 


10/74 


OOZ 21 




1/96 


5/74 


4/51 


4/45 


7/50 


12Z 22 




4/87 


5/61 


5/48 


7/42 


6/50 


12Z 23 




0/90 


5/76 


-6/73 


-3/71 


5/68 


OOZ 24 




-5/95 


0/65 


6/63 


5/69 


0/74 


12Z 24 




-3/97 


-4/77 


-7/71 


-2/70 


4/66 


Table II 


Pillow and RMSE 


of persistence verifications 


at 5 level 




for 


a 24-hr forecast during November 


1962. 




Initial 


Time 






Level 










300mb 


500mb 


700mb 


850mb 


lOOOmb 


12Z 19 




-1/143 


-3/96 


2/77 


0/89 


-11/106 


OOZ 20 




-35/122 


-13/95 


0/92 


2/106 


-1/116 


12Z 21 




-16/135 


1/100 


7/90 


12/94 


18/89 


12Z 22 




0/122 


10/91 


14/82 


13/83 


14/91 
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Table III. Pillow and RMSE of persistence verifications at 5 levels 
for 36, 48, 60 and 72-hr forecasts during November 1962, 



Initial Forecast 



Time 


Interval 






Level 










300mb 


500rab 


700mb 


850mb 


lOOOmb 


12Z 19 


36 


-13/159 


-6/115 


0/116 


0/130 


0/155 


12Z 22 


36 


0/139 


5/140 


7/134 


10/133 


20/139 


12Z 19 


48 


-11/188 


-1/147 


4/137 


4/150 


6/171 


12Z 22 


48 


-4/144 


6/151 


14/163 


16/174 


20/185 


12Z 22 


60 


-8/178 


1/188 


7/203 


13/208 


24/213 


12Z 22 


72 


-17/231 


4/223 


20/223 


21/227 


23/233 
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